##########################################################
# Código para os modelos logísticos por tipo institucional
# 
# Autoria: Leonardo Rodrigues
# Código escrito e executado no RStudio (Version 1.3.1073)
# Tempo aproximado para rodar o script: 2.75 horas
##########################################################
start_time <- Sys.time()

### 1 Preparação ----
# Leitura dos pacotes

library(webshot)
library(dplyr)
library(stringr)
library(tidyr)
library(tools)
library(summarytools)
library(stargazer)
library(mice)
library(miceadds)
library(nnet)
library(modelsummary)

#Leitura do banco
enade <- read.table("data/originaisEnade/enade.txt", sep="\t")

#Ajustes das variáveis
table(enade$cor)

#Mudando nome das categorias para aproveitar script anterior
enade$cor = ifelse(enade$cor == "Brancos", "Branca",
                   ifelse(enade$cor == "Negros", "Não-Branca", NA))

#renomeando variáveis para aproveitar script anterior
#estabelecendo a categoria de referência para os modelos
enade$year <- as.factor(enade$ciclo)
enade$year = relevel(enade$year, ref = "2")
enade$tese = ifelse(str_detect(enade$curso,"ENGENHARIA"), enade$curso, "outros")
enade11 <- enade
enade11$year <- enade11$ciclo
enade11$year <- as.factor(enade11$year)
enade11$year = relevel(enade11$year, ref = "2")
enade11$grupo2 <- enade11$diploma
enade11$grupo2 <- as.factor(enade11$grupo2)
enade11$grupo2 = relevel(enade11$grupo2, ref = "Engenharias")
enade11$regiao <- enade11$regiao1
rm(enade)

#selecionando as variáveis de interesse
enade11 <- select(enade11, grupo2, age, regiao, sexo, cor, edu, year, inst)

#Estabelecendo as categorias de referência para cada variável
enade11$grupo2 = relevel(enade11$grupo2, ref = "Engenharias")
enade11$age <- as.factor(enade11$age)
enade11$age = relevel(enade11$age, ref = "20-30")
enade11$regiao <- as.factor(enade11$regiao)
enade11$regiao = relevel(enade11$regiao, ref = "Sul e Sudeste")
enade11$sexo <- as.factor(enade11$sexo)
enade11$sexo = relevel(enade11$sexo, ref = "Masculino")
enade11$cor <- as.factor(enade11$cor)
enade11$cor = relevel(enade11$cor, ref = "Branca")
enade11$edu <- as.factor(enade11$edu)
enade11$edu = relevel(enade11$edu, ref = "Ensino superior ou mais")
enade11$inst <- as.factor(enade11$inst)
enade11$inst = relevel(enade11$inst, ref = "Pública")

### 2 Imputação de valores para os missing ----
#Separando os bancos por ciclo
enade1 <- enade11[enade11$year==2,]
enade2 <- enade11[enade11$year==4,]
rm(enade11)

#Imputação para cada ciclo
imputed_data1 <- mice(enade1, m=5, maxit=25, method = 'logreg')

imputed_data2 <- mice(enade2, m=5, maxit=25, method = 'logreg')

#Junção dos ciclos
imputed_total <- rbind(imputed_data1, imputed_data2) #juntando os bancos de dados imputados
rm(imputed_data1, imputed_data2)

### 3 Modelos logísticos para cada área por tipo institucional ----

#Modelo para as engenharias
#separar os bancos imputados por área de formação
enade <- mice::complete(imputed_total,"long",include = T)

eng <- enade[which(enade$grupo2 == "Engenharias"), ]
eng <- as.mids(eng)

#modelo
model <- with(eng, glm(inst ~ age + regiao + sexo*year + cor*year + edu*year, family=binomial(link='logit')))
rm(eng)

model1 <- pool(model) #juntando os cinco modelos

#substituir "estimate" pelo "Pooled Estimate"
coeficientes <- model$analyses[[1]]
coeficientes$coefficients[[1]] <- summary(model1)$estimate[1]
coeficientes$coefficients[[2]] <- summary(model1)$estimate[2]
coeficientes$coefficients[[3]] <- summary(model1)$estimate[3]
coeficientes$coefficients[[4]] <- summary(model1)$estimate[4]
coeficientes$coefficients[[5]] <- summary(model1)$estimate[5]
coeficientes$coefficients[[6]] <- summary(model1)$estimate[6]
coeficientes$coefficients[[7]] <- summary(model1)$estimate[7]
coeficientes$coefficients[[8]] <- summary(model1)$estimate[8]
coeficientes$coefficients[[9]] <- summary(model1)$estimate[9]
coeficientes$coefficients[[10]] <- summary(model1)$estimate[10]

#Salvando o resumo do pooled model
summaryEng <- msummary(model1, stars = TRUE)
readr::write_file(summaryEng, 'data/output/summaryEng.html')

#obtendo a razão entre as probabilidades preditas com a estratégia semelhante ao modelo anterior
pps <- data.frame("area" = c("Engenharia", "Medicina","Direito",
                             "Bacharelados","Licenciatura","Tecnólogos"),
                  "MF1" = NA, "SM1" = NA, "BN1" = NA,
                  "MF2" = NA, "SM2" = NA, "BN2" = NA)

#educação dos pais
mediaEDU <- data.frame(edu = c("Ensino superior ou mais","Até Ensino Médio", "Ensino superior ou mais", "Até Ensino Médio"),
                       year = c("2", "2", "4", "4"),
                       sexo = c("Masculino", "Masculino", "Masculino", "Masculino"),
                       cor = c("Branca","Branca", "Branca","Branca"),
                       age = c("20-30", "20-30", "20-30", "20-30"),
                       regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaEDU[, c("pred.prob")] <- predict(coeficientes, newdata=mediaEDU, type="response")



#sexo

mediaSEXO <- data.frame(edu = c("Até Ensino Médio","Até Ensino Médio", "Até Ensino Médio", "Até Ensino Médio"),
                        year = c("2", "2", "4", "4"),
                        sexo = c("Masculino", "Feminino", "Masculino", "Feminino"),
                        cor = c("Branca","Branca", "Branca","Branca"),
                        age = c("20-30", "20-30", "20-30", "20-30"),
                        regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaSEXO[, c("pred.prob")] <- predict(coeficientes, newdata=mediaSEXO, type="response")


#cor/raça

mediaCOR <- data.frame(edu = c("Até Ensino Médio","Até Ensino Médio", "Até Ensino Médio", "Até Ensino Médio"),
                       year = c("2", "2", "4", "4"),
                       sexo = c("Masculino", "Masculino", "Masculino", "Masculino"),
                       cor = c("Branca","Não-Branca", "Branca","Não-Branca"),
                       age = c("20-30", "20-30", "20-30", "20-30"),
                       regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaCOR[, c("pred.prob")] <- predict(coeficientes, newdata=mediaCOR, type="response")


pps[pps$area == "Engenharia", "SM1"] <- mediaEDU[1,7]/mediaEDU[2,7]
pps[pps$area == "Engenharia", "SM2"] <- mediaEDU[3,7]/mediaEDU[4,7]
pps[pps$area == "Engenharia", "MF1"] <- mediaSEXO[1,7]/mediaSEXO[2,7]
pps[pps$area == "Engenharia", "MF2"] <- mediaSEXO[3,7]/mediaSEXO[4,7]
pps[pps$area == "Engenharia", "BN1"] <- mediaCOR[1,7]/mediaCOR[2,7]
pps[pps$area == "Engenharia", "BN2"] <- mediaCOR[3,7]/mediaCOR[4,7]

#Modelo para a medicina
med <- enade[which(enade$grupo2 == "Medicina"), ]

med <- as.mids(med)

#modelo
model <- with(med, glm(inst ~ age + regiao + sexo*year + cor*year + edu*year, family=binomial(link='logit')))
rm(med)

model1 <- pool(model) #juntando os cinco modelos

#substituir "estimate" pelo "Pooled Estimate"
coeficientes <- model$analyses[[1]]
coeficientes$coefficients[[1]] <- summary(model1)$estimate[1]
coeficientes$coefficients[[2]] <- summary(model1)$estimate[2]
coeficientes$coefficients[[3]] <- summary(model1)$estimate[3]
coeficientes$coefficients[[4]] <- summary(model1)$estimate[4]
coeficientes$coefficients[[5]] <- summary(model1)$estimate[5]
coeficientes$coefficients[[6]] <- summary(model1)$estimate[6]
coeficientes$coefficients[[7]] <- summary(model1)$estimate[7]
coeficientes$coefficients[[8]] <- summary(model1)$estimate[8]
coeficientes$coefficients[[9]] <- summary(model1)$estimate[9]
coeficientes$coefficients[[10]] <- summary(model1)$estimate[10]

#Salvando o resumo do pooled model
summaryMed <- msummary(model1, stars = TRUE)
readr::write_file(summaryMed, 'data/output/summaryMed.html')

#obtendo a razão entre as probabilidades preditas com a estratégia semelhante ao modelo anterior
#educação dos pais
mediaEDU <- data.frame(edu = c("Ensino superior ou mais","Até Ensino Médio", "Ensino superior ou mais", "Até Ensino Médio"),
                       year = c("2", "2", "4", "4"),
                       sexo = c("Feminino", "Feminino", "Feminino", "Feminino"),
                       cor = c("Branca","Branca", "Branca","Branca"),
                       age = c("20-30", "20-30", "20-30", "20-30"),
                       regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaEDU[, c("pred.prob")] <- predict(coeficientes, newdata=mediaEDU, type="response")

#sexo

mediaSEXO <- data.frame(edu = c("Ensino superior ou mais","Ensino superior ou mais", "Ensino superior ou mais", "Ensino superior ou mais"),
                        year = c("2", "2", "4", "4"),
                        sexo = c("Masculino", "Feminino", "Masculino", "Feminino"),
                        cor = c("Branca","Branca", "Branca","Branca"),
                        age = c("20-30", "20-30", "20-30", "20-30"),
                        regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaSEXO[, c("pred.prob")] <- predict(coeficientes, newdata=mediaSEXO, type="response")

#cor

mediaCOR <- data.frame(edu = c("Ensino superior ou mais","Ensino superior ou mais", "Ensino superior ou mais", "Ensino superior ou mais"),
                       year = c("2", "2", "4", "4"),
                       sexo = c("Masculino", "Masculino", "Masculino", "Masculino"),
                       cor = c("Branca","Não-Branca", "Branca","Não-Branca"),
                       age = c("20-30", "20-30", "20-30", "20-30"),
                       regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaCOR[, c("pred.prob")] <- predict(coeficientes, newdata=mediaCOR, type="response")

pps[pps$area == "Medicina", "SM1"] <- mediaEDU[1,7]/mediaEDU[2,7]
pps[pps$area == "Medicina", "SM2"] <- mediaEDU[3,7]/mediaEDU[4,7]
pps[pps$area == "Medicina", "MF1"] <- mediaSEXO[1,7]/mediaSEXO[2,7]
pps[pps$area == "Medicina", "MF2"] <- mediaSEXO[3,7]/mediaSEXO[4,7]
pps[pps$area == "Medicina", "BN1"] <- mediaCOR[1,7]/mediaCOR[2,7]
pps[pps$area == "Medicina", "BN2"] <- mediaCOR[3,7]/mediaCOR[4,7]

#Modelo para o direito
dir <- enade[which(enade$grupo2 == "Direito"), ]

dir <- as.mids(dir)

#modelo
model <- with(dir, glm(inst ~ age + regiao + sexo*year + cor*year + edu*year, family=binomial(link='logit')))
rm(dir)

model1 <- pool(model) #juntando os cinco modelos

#substituir "estimate" pelo "Pooled Estimate"
coeficientes <- model$analyses[[1]]
coeficientes$coefficients
coeficientes$coefficients[[1]] <- summary(model1)$estimate[1]
coeficientes$coefficients[[2]] <- summary(model1)$estimate[2]
coeficientes$coefficients[[3]] <- summary(model1)$estimate[3]
coeficientes$coefficients[[4]] <- summary(model1)$estimate[4]
coeficientes$coefficients[[5]] <- summary(model1)$estimate[5]
coeficientes$coefficients[[6]] <- summary(model1)$estimate[6]
coeficientes$coefficients[[7]] <- summary(model1)$estimate[7]
coeficientes$coefficients[[8]] <- summary(model1)$estimate[8]
coeficientes$coefficients[[9]] <- summary(model1)$estimate[9]
coeficientes$coefficients[[10]] <- summary(model1)$estimate[10]

#Salvando o resumo do pooled model
summaryDir <- msummary(model1, stars = TRUE)
readr::write_file(summaryDir, 'data/output/summaryDir.html')

#obtendo a razão entre as probabilidades preditas com a estratégia semelhante ao modelo anterior
#educação dos pais
mediaEDU <- data.frame(edu = c("Ensino superior ou mais","Até Ensino Médio", "Ensino superior ou mais", "Até Ensino Médio"),
                       year = c("2", "2", "4", "4"),
                       sexo = c("Feminino", "Feminino", "Feminino", "Feminino"),
                       cor = c("Branca","Branca", "Branca","Branca"),
                       age = c("20-30", "20-30", "20-30", "20-30"),
                       regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaEDU[, c("pred.prob")] <- predict(coeficientes, newdata=mediaEDU, type="response")

#sexo

mediaSEXO <- data.frame(edu = c("Ensino superior ou mais","Ensino superior ou mais", "Ensino superior ou mais", "Ensino superior ou mais"),
                        year = c("2", "2", "4", "4"),
                        sexo = c("Masculino", "Feminino", "Masculino", "Feminino"),
                        cor = c("Branca","Branca", "Branca","Branca"),
                        age = c("20-30", "20-30", "20-30", "20-30"),
                        regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaSEXO[, c("pred.prob")] <- predict(coeficientes, newdata=mediaSEXO, type="response")

#cor

mediaCOR <- data.frame(edu = c("Ensino superior ou mais","Ensino superior ou mais", "Ensino superior ou mais", "Ensino superior ou mais"),
                       year = c("2", "2", "4", "4"),
                       sexo = c("Feminino", "Feminino", "Feminino", "Feminino"),
                       cor = c("Branca","Não-Branca", "Branca","Não-Branca"),
                       age = c("20-30", "20-30", "20-30", "20-30"),
                       regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaCOR[, c("pred.prob")] <- predict(coeficientes, newdata=mediaCOR, type="response")

pps[pps$area == "Direito", "SM1"] <- mediaEDU[1,7]/mediaEDU[2,7]
pps[pps$area == "Direito", "SM2"] <- mediaEDU[3,7]/mediaEDU[4,7]
pps[pps$area == "Direito", "MF1"] <- mediaSEXO[1,7]/mediaSEXO[2,7]
pps[pps$area == "Direito", "MF2"] <- mediaSEXO[3,7]/mediaSEXO[4,7]
pps[pps$area == "Direito", "BN1"] <- mediaCOR[1,7]/mediaCOR[2,7]
pps[pps$area == "Direito", "BN2"] <- mediaCOR[3,7]/mediaCOR[4,7]

#Modelo para os bacharelados

bac <- enade[which(enade$grupo2 == "Bacharelados"), ]

bac <- as.mids(bac)

#modelo
model <- with(bac, glm(inst ~ age + regiao + sexo*year + cor*year + edu*year, family=binomial(link='logit')))
rm(bac)

model1 <- pool(model) #juntando os cinco modelos

#substituir "estimate" pelo "Pooled Estimate"
coeficientes <- model$analyses[[1]]
coeficientes$coefficients
coeficientes$coefficients[[1]] <- summary(model1)$estimate[1]
coeficientes$coefficients[[2]] <- summary(model1)$estimate[2]
coeficientes$coefficients[[3]] <- summary(model1)$estimate[3]
coeficientes$coefficients[[4]] <- summary(model1)$estimate[4]
coeficientes$coefficients[[5]] <- summary(model1)$estimate[5]
coeficientes$coefficients[[6]] <- summary(model1)$estimate[6]
coeficientes$coefficients[[7]] <- summary(model1)$estimate[7]
coeficientes$coefficients[[8]] <- summary(model1)$estimate[8]
coeficientes$coefficients[[9]] <- summary(model1)$estimate[9]
coeficientes$coefficients[[10]] <- summary(model1)$estimate[10]

#Salvando o resumo do pooled model
summaryBac <- msummary(model1, stars = TRUE)
readr::write_file(summaryBac, 'data/output/summaryBac.html')

#obtendo a razão entre as probabilidades preditas com a estratégia semelhante ao modelo anterior
#educação dos pais
mediaEDU <- data.frame(edu = c("Ensino superior ou mais","Até Ensino Médio", "Ensino superior ou mais", "Até Ensino Médio"),
                       year = c("2", "2", "4", "4"),
                       sexo = c("Feminino", "Feminino", "Feminino", "Feminino"),
                       cor = c("Branca","Branca", "Branca","Branca"),
                       age = c("20-30", "20-30", "20-30", "20-30"),
                       regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaEDU[, c("pred.prob")] <- predict(coeficientes, newdata=mediaEDU, type="response")


#sexo

mediaSEXO <- data.frame(edu = c("Até Ensino Médio","Até Ensino Médio", "Até Ensino Médio", "Até Ensino Médio"),
                        year = c("2", "2", "4", "4"),
                        sexo = c("Masculino", "Feminino", "Masculino", "Feminino"),
                        cor = c("Branca","Branca", "Branca","Branca"),
                        age = c("20-30", "20-30", "20-30", "20-30"),
                        regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaSEXO[, c("pred.prob")] <- predict(coeficientes, newdata=mediaSEXO, type="response")


#cor

mediaCOR <- data.frame(edu = c("Até Ensino Médio","Até Ensino Médio", "Até Ensino Médio", "Até Ensino Médio"),
                       year = c("2", "2", "4", "4"),
                       sexo = c("Feminino", "Feminino", "Feminino", "Feminino"),
                       cor = c("Branca","Não-Branca", "Branca","Não-Branca"),
                       age = c("20-30", "20-30", "20-30", "20-30"),
                       regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaCOR[, c("pred.prob")] <- predict(coeficientes, newdata=mediaCOR, type="response")

pps[pps$area == "Bacharelados", "SM1"] <- mediaEDU[1,7]/mediaEDU[2,7]
pps[pps$area == "Bacharelados", "SM2"] <- mediaEDU[3,7]/mediaEDU[4,7]
pps[pps$area == "Bacharelados", "MF1"] <- mediaSEXO[1,7]/mediaSEXO[2,7]
pps[pps$area == "Bacharelados", "MF2"] <- mediaSEXO[3,7]/mediaSEXO[4,7]
pps[pps$area == "Bacharelados", "BN1"] <- mediaCOR[1,7]/mediaCOR[2,7]
pps[pps$area == "Bacharelados", "BN2"] <- mediaCOR[3,7]/mediaCOR[4,7]

#licenciatura

lic <- enade[which(enade$grupo2 == "Licenciaturas"), ]

lic <- as.mids(lic)

#modelo
model <- with(lic, glm(inst ~ age + regiao + sexo*year + cor*year + edu*year, family=binomial(link='logit')))
rm(lic)

model1 <- pool(model) #juntando os cinco modelos

#substituir "estimate" pelo "Pooled Estimate"
coeficientes <- model$analyses[[1]]
coeficientes$coefficients[[1]] <- summary(model1)$estimate[1]
coeficientes$coefficients[[2]] <- summary(model1)$estimate[2]
coeficientes$coefficients[[3]] <- summary(model1)$estimate[3]
coeficientes$coefficients[[4]] <- summary(model1)$estimate[4]
coeficientes$coefficients[[5]] <- summary(model1)$estimate[5]
coeficientes$coefficients[[6]] <- summary(model1)$estimate[6]
coeficientes$coefficients[[7]] <- summary(model1)$estimate[7]
coeficientes$coefficients[[8]] <- summary(model1)$estimate[8]
coeficientes$coefficients[[9]] <- summary(model1)$estimate[9]
coeficientes$coefficients[[10]] <- summary(model1)$estimate[10]

#Salvando o resumo do pooled model
summaryLic <- msummary(model1, stars = TRUE)
readr::write_file(summaryLic, 'data/output/summaryLic.html')

#obtendo a razão entre as probabilidades preditas com a estratégia semelhante ao modelo anterior
#educação dos pais

mediaEDU <- data.frame(edu = c("Ensino superior ou mais","Até Ensino Médio", "Ensino superior ou mais", "Até Ensino Médio"),
                       year = c("2", "2", "4", "4"),
                       sexo = c("Feminino", "Feminino", "Feminino", "Feminino"),
                       cor = c("Não-Branca","Não-Branca", "Não-Branca","Não-Branca"),
                       age = c("20-30", "20-30", "20-30", "20-30"),
                       regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaEDU[, c("pred.prob")] <- predict(coeficientes, newdata=mediaEDU, type="response")


#sexo

mediaSEXO <- data.frame(edu = c("Até Ensino Médio","Até Ensino Médio", "Até Ensino Médio", "Até Ensino Médio"),
                        year = c("2", "2", "4", "4"),
                        sexo = c("Masculino", "Feminino", "Masculino", "Feminino"),
                        cor = c("Não-Branca","Não-Branca", "Não-Branca","Não-Branca"),
                        age = c("20-30", "20-30", "20-30", "20-30"),
                        regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaSEXO[, c("pred.prob")] <- predict(coeficientes, newdata=mediaSEXO, type="response")

#cor

mediaCOR <- data.frame(edu = c("Até Ensino Médio","Até Ensino Médio", "Até Ensino Médio", "Até Ensino Médio"),
                       year = c("2", "2", "4", "4"),
                       sexo = c("Feminino", "Feminino", "Feminino", "Feminino"),
                       cor = c("Branca","Não-Branca", "Branca","Não-Branca"),
                       age = c("20-30", "20-30", "20-30", "20-30"),
                       regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaCOR[, c("pred.prob")] <- predict(coeficientes, newdata=mediaCOR, type="response")

pps[pps$area == "Licenciatura", "SM1"] <- mediaEDU[1,7]/mediaEDU[2,7]
pps[pps$area == "Licenciatura", "SM2"] <- mediaEDU[3,7]/mediaEDU[4,7]
pps[pps$area == "Licenciatura", "MF1"] <- mediaSEXO[1,7]/mediaSEXO[2,7]
pps[pps$area == "Licenciatura", "MF2"] <- mediaSEXO[3,7]/mediaSEXO[4,7]
pps[pps$area == "Licenciatura", "BN1"] <- mediaCOR[1,7]/mediaCOR[2,7]
pps[pps$area == "Licenciatura", "BN2"] <- mediaCOR[3,7]/mediaCOR[4,7]

#obtendo a razão entre as probabilidades preditas com a estratégia semelhante ao modelo anterior
#tecnológicos
tec <- enade[which(enade$grupo2 == "Tecnólogos"), ]

tec <- as.mids(tec)

#modelo
model <- with(tec, glm(inst ~ age + regiao + sexo*year + cor*year + edu*year, family=binomial(link='logit')))
rm(tec)

model1 <- pool(model) #juntando os cinco modelos

#substituir "estimate" pelo "Pooled Estimate"
coeficientes <- model$analyses[[1]]
coeficientes$coefficients[[1]] <- summary(model1)$estimate[1]
coeficientes$coefficients[[2]] <- summary(model1)$estimate[2]
coeficientes$coefficients[[3]] <- summary(model1)$estimate[3]
coeficientes$coefficients[[4]] <- summary(model1)$estimate[4]
coeficientes$coefficients[[5]] <- summary(model1)$estimate[5]
coeficientes$coefficients[[6]] <- summary(model1)$estimate[6]
coeficientes$coefficients[[7]] <- summary(model1)$estimate[7]
coeficientes$coefficients[[8]] <- summary(model1)$estimate[8]
coeficientes$coefficients[[9]] <- summary(model1)$estimate[9]
coeficientes$coefficients[[10]] <- summary(model1)$estimate[10]

#Salvando o resumo do pooled model
summaryTec <- msummary(model1, stars = TRUE)
readr::write_file(summaryTec, 'data/output/summaryTec.html')

#educação dos pais
mediaEDU <- data.frame(edu = c("Ensino superior ou mais","Até Ensino Médio", "Ensino superior ou mais", "Até Ensino Médio"),
                       year = c("2", "2", "4", "4"),
                       sexo = c("Feminino", "Feminino", "Feminino", "Feminino"),
                       cor = c("Branca","Branca", "Branca","Branca"),
                       age = c("20-30", "20-30", "20-30", "20-30"),
                       regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaEDU[, c("pred.prob")] <- predict(coeficientes, newdata=mediaEDU, type="response")


#sexo

mediaSEXO <- data.frame(edu = c("Até Ensino Médio","Até Ensino Médio", "Até Ensino Médio", "Até Ensino Médio"),
                        year = c("2", "2", "4", "4"),
                        sexo = c("Masculino", "Feminino", "Masculino", "Feminino"),
                        cor = c("Branca","Branca", "Branca","Branca"),
                        age = c("20-30", "20-30", "20-30", "20-30"),
                        regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaSEXO[, c("pred.prob")] <- predict(coeficientes, newdata=mediaSEXO, type="response")


#cor

mediaCOR <- data.frame(edu = c("Até Ensino Médio","Até Ensino Médio", "Até Ensino Médio", "Até Ensino Médio"),
                       year = c("2", "2", "4", "4"),
                       sexo = c("Feminino", "Feminino", "Feminino", "Feminino"),
                       cor = c("Branca","Não-Branca", "Branca","Não-Branca"),
                       age = c("20-30", "20-30", "20-30", "20-30"),
                       regiao = c("Sul e Sudeste","Sul e Sudeste", "Sul e Sudeste","Sul e Sudeste"))

mediaCOR[, c("pred.prob")] <- predict(coeficientes, newdata=mediaCOR, type="response")

pps[pps$area == "Tecnólogos", "SM1"] <- mediaEDU[1,7]/mediaEDU[2,7]
pps[pps$area == "Tecnólogos", "SM2"] <- mediaEDU[3,7]/mediaEDU[4,7]
pps[pps$area == "Tecnólogos", "MF1"] <- mediaSEXO[1,7]/mediaSEXO[2,7]
pps[pps$area == "Tecnólogos", "MF2"] <- mediaSEXO[3,7]/mediaSEXO[4,7]
pps[pps$area == "Tecnólogos", "BN1"] <- mediaCOR[1,7]/mediaCOR[2,7]
pps[pps$area == "Tecnólogos", "BN2"] <- mediaCOR[3,7]/mediaCOR[4,7]


#salvando os resultados para fazer o gráfico
write.table(pps, "data/output/ppslogit.txt", sep="\t")

### END OF SCRIPT ###
end_time <- Sys.time()
end_time - start_time #2.750452 